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Inhomogeneous  conditions  at  open  boundaries  for  wave  propagation  problems 


1.  Introduction 

When  wave  propagation  problems  are  solved  numerically,  it  is  often  necessary  to  introduce  artificial 
boundaries.  In  most  cases  no  data  are  available  at  these  boundaries,  and  therefore  it  is  necessary  to 
construct  boundary  conditions  which  in  some  way  accounts  for  the  behavior  of  the  solution  outside 
the  computational  domain  -fit  Hence  some  assumption  must  be  made  which  makes  it  possible  to 
solve  the  problem  exactly  outside  or  alternatively  to  compute  the  solution  approximately  in  a 
simple  way.  One  such  assumption  is  that  the  waves  are  propagating  only  in  the  outward  direction 
across  the  boundary.  This  is  the  basis  for  various  procedures  in  common  use,  the  most  general 
class  being  the  absorbing  boundary  condition^  constructed  by  Engquist  and  Majda  [1],  [2],  see 
also  Higdon  [7].  In  [4]  and  [5]  more  general  procedures  were  discussed,  in  particular  the  case  with 
non-zero  initial  data  outside  the  computational  domain  was  considered. 

^  The  absorbing  boundary  conditions  of  higher  order  are  formulated  in  terms  of  differentiated 
function^ As  indicated  in  [4],^ this  leads  necessarily  to  weakly  ill-posed  problems,  and  as  a  conse¬ 
quence  also  to  unstable  numerical  methods.  This  has  been  pointed  out  and  discussed  by  Higdon  {7}  * 
for  the  scalar  wave-equatiodT The  author  argues  that  for  conditions  of  order  two  or  less,  the  possible 
effects  of  the  instability  is  outweighed  by  the  small  reflection  coefficients. 

Tn  this  paper  we  sha^  investigate'^ his  problem  further.  When  the  initial  data  do  not  have 
compact  support  within  the  computational  domain,  or  when  some  outside  source  is  present,  the 
boundary  conditions  become  inhomogeneous.  This  tends  to  amplify  the  influence  of  the  instability, 
and  it  becomes  essential  to  modify  the  boundary  procedure..  We  shall  demonstrate  how  this  can 
be  done. 

In  Section  2  we  treat  the  scalar  wave  equation,  and  we  demonstrate  the  effect  of  higher  order 
derivatives  in  the  boundary  conditions.  More  smoothness  of  the  solution  is  required  to  compensate 
for  the  weak  ill-posedness,  and  the  approximation  becomes  more  sensitive  to  perturbations.  We 
also  derive  a  stability  estimate  for  the  usual  five-point  approximation  in  one  space  dimension  with 
the  inhomogeneous  version  of  the  first  absorbing  boundary  condition  suggested  in  [1].  The  estimate 
is  an  exact  discrete  version  of  the  one  obtained  for  the  continuous  problem. 

In  Section  3  we  discuss  general  first  order  hyperbolic  systems,  and  the  implementation  of 
higher  order  accurate  boundary  conditions.  These  are  constructed  such  that  the  weak  ill-posedness 
is  removed,  which  makes  it  possible  to  construct  stable  approximations. 

In  a  recent  paper,  Howell  and  Trefethen  [8]  have  shown  ill-posedness  for  migration  equations 
in  connection  with  certain  boundary  conditions.  In  this  paper  we  consider  the  wave  propagation 
problems  in  its  original  form,  i.e.,  the  differential  equations  are  hyperbolic. 

No  attempt  is  made  here  to  discuss  various  considerations  concerning  reflection  properties.  This 
has  been  done  in  the  papers  mentioned  above.  The  purpose  of  this  paper  is  only  to  illuminate  the 
influence  of  the  weak  instabilities  inherent  with  higher  order  derivatives  in  the  boundary  conditions, 
and  how  to  remove  them.  These  principles  can  also  be  applied  to  problems  which  are  not  of  pure 
wave-  propagation  type. 
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2.  The  wave  equation 

We  consider  the  wave  equation  where  a  scaling  has  been  done  such  that  the  propagation  speed 
is  one.  For  convenience,  we  treat  only  the  quarter-space  problem 

{fat  =  <t>*x  +  4>w  ,  0  <  x  <  oo,  -oo  <  y  <  oo,  0  <  t, 

fax,Q)  =  f(x,y), 
fa(x,0)  =  d(x,y), 

where  we  need  also  a  boundary  condition  at  x  =  0.  The  general  theory  for  hyperbolic  initial¬ 
boundary  value  problems  is  given  in  [9].  It  is  formulated  for  first  order  systems,  but  the  same 
technique  (the  normal  mode  analysis)  for  analyzing  well-posedness  can  be  applied  to  differential 
equations  and  boundary  conditions  where  higher  order  derivatives  occur.  We  remark  further  on 
this  in  Section  3. 

The  problem  is  Laplace- transformed  in  the  f-direction  and  Fourier- transformed  in  the  y- 
direction,  s  and  u  being  the  dual  variables  respectively.  For  /  =  d  =  0,  the  transformed  solution 
to  (2.1)  is 

fax,u,s)  =  e-‘'/^faO,uJs), 

where  the  square  root  is  defined  such  that  Re  \Jz  >  0  for  all  z.  If  the  transformed  boundary 
condition  is 

fl(u,s)<£(0,u;,3)  =  j(w,s) , 

the  uniform  Kreiss  condition  required  for  well-posedness  is 

(2.2)  |f?(w,  s)|  >  6  >  0 ,  Re  s  >  0,  Im  w  =  0 . 

When  the  wave  function  <t>  is  specified  at  *  =  0,  such  that 

(2.3)  faO,y,t)  =  g(y,t), 

we  have  B  =  1,  and  the  condition  (2.2)  is  trivially  satisfied. 

Consider  next  the  boundary  condition 

(2.4)  fa  -  <t>*  =  g(y,t) ,  x  =  0, 

which  is  the  inhomogeneous  version  of  the  first  absorbing  boundary  condition.  The  problem  of 
providing  the  data  g(y,  t)  will  be  discussed  in  Section  3.  For  (2.4)  we  have 

B(u,a)  =  3+  y/s1  +  , 

and  obviously  the  condition  (2.2)  is  violated  foru;  =  s  =  0.  When  the  Kreiss  condition  is  violated  for 
a  =  so  on  the  imaginary  axis,  so  is  called  a  generalized  eigenvalue.  In  such  a  case  the  problem  is  said 
to  be  weakly  ill-posed.  However,  it  seems  natural  to  accept  first  order  derivatives  in  the  boundary 
conditions  when  the  differential  equation  is  of  second  order,  and  we  shall  further  investigate  this 


Since  the  generalized  eigenvalue  occurs  for  u>  =  0,  it  is  sufficient  to  study  the  one-dimensional 
case.  We  consider  the  problem 

{(a)  4>tt  =  <j>xx,  0  <  x  <  oo,  0  <  t 

(b)  4>t  -  4>x  =  g(t) ,  *  =  0 

(c)  *  =  /(z),  t  =  0 

(d)  <j>t  =  d(x),  t  =  0. 

It  is  assumed  that  /,  d  have  compact  support,  such  that  for  any  given  t,  4>{x,  t)  also  has  compact 
support.  The  solution  of  (2.5)  can  be  given  explicitly,  but  we  derive  the  general  estimate  for 
the  solution  in  order  to  get  a  procedure  that  can  be  modified  in  a  straightforward  manner  to 
discretizations  of  (2.5). 

The  wave  equation  is  rewritten  as  a  first  order  system.  With 

«  =  <t>t ,  v  =  4>x 

we  get 


rui  TO  ll  rui 

UJi  “  [l  OJ  LJ,’ 

0<z<oo,  0<f 

ll 

& 

i 

3 

z  =  0 

u  =  d(z) , 

t  =  0 

V  =  /'(z), 

t  =  0. 

Since  the  ingoing  characteristic  variable  is  prescribed,  we  get  aw  estimate  also  for  the  boundary 
values  of  the  solution.  A  straightforward  application  of  the  energy  method  gives 

MOII2  +  Ht)\\7  +  /‘(|«(0,r)|J  +  |r(0,r)|’)dr 
(2.6)  t 

<C'(||d||I  +  ||ri|J  +  j[  Is(i-)!2*-) 

where  ||  •  ||  denotes  the  L2-norm.  Note  that  the  reformulation  as  a  first  order  system  has  eliminated 
the  generalized  eigenvalue  a  =  0,  making  the  estimate  (2.6)  possible. 

In  order  to  get  an  estimate  also  for  0,  we  use  the  fact  that  <p{x,t)  has  compact  support,  and 
obtain  for  some  zo  =  xq (t), 

mm2  =  j‘°  dx  <  c(iiv(t)ii2  +  i<£(o,t)i2) , 

(2.7)  WO, 0I2  =  1/(0)  +  J*  u(0,  r)  dr|J  <  2  (|/(0)|2  +  J‘  |u(0,  r)|2  dr)  . 

The  integral  in  the  right  hand  side  of  (2.7)  can  be  estimated  by  using  (2.6),  and  we  arrive  at  the 
final  estimate 

IW0II*  +  IWWII*  +  IW(0l|2  +  W0, t)|2  +  /‘(|^(0,r)|2  +  W(0,r)|2)  dr 

Jo 


<  C  (||d||2  +  ll/'H2  +  |/(0)|2  +  jf  |ff(r)|2  dr)  , 
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where  C  depends  on  t.  Here  the  derivatives  of  0  and  of  the  initial  data  are  included,  which  is 
very  natural.  These  derivatives  are  physical  quantities  of  prime  interest,  and  (2.8)  is  therefore  a 
useful  estimate.  For  non-smooth  initial  data  an  estimate  can  still  be  obtained.  With  the  variable 
substitution 

it  =  <(> 


=  J*Mx,t)  dr  +  J’  d(t)dt 


we  get  after  integrating  (2.5b)  with  respect  to  t, 


a-p  i]  ci.. « 


i  <  oo,  0  <  t 


u-v  =  g(t) , 

«  =  /(*) , 
v  =  tid(Z)dt, 


x  =  0 


t  =  0, 


where 


§(t)  =  /  9(r)dr  +  f( 0) . 
Jo 


Corresponding  to  (2.6)  we  get 


IWOII 


1  +  J*  |«0,  r)!1  dr  <C  (||/||J  +  ||  J‘  d( O  d*||J  +  J'\9{T)\2  dr) 


Next,  consider  the  boundary  condition 


<l>it  -  4*t  ~  =  $(iM)- 


For  g(y,  t)  =  0,  this  is  the  second  absorbing  boundary  condition.  In  transformed  space,  the  Kreiss 
condition  obtained  from  (2.9)  is 


(2.10) 


B(w,i)  s  «(i  +  Vj’Tw^)  +  /  0 ,  Rea>0. 


For  ur  =  0,  a  =  0,  the  condition  is  not  fulfilled,  and  we  again  have  a  generalized  eigenvalue.  (This 
case  was  not  included  in  the  analysis  carried  out  in  [1],  [2],  and  [10].) 

In  the  analysis  of  the  one-dimensional  problem  above,  the  final  estimate  was  obtained  by 
considering  a  first  order  system  for  the  derivatives  of  <f>.  For  the  higher  order  boundary  condition, 
we  apply  the  same  technique  once  more.  The  functions  u  =  4>u,  v  =  4>tt  satisfy  the  same  differential 
equation  as  u,  v  above,  but  with  new  initial  conditions,  which  can  be  derived  by  using  the  differential 
equation.  We  get 


«]  _  ro  i]  [«■ 

J  * _  [i  UJj 


,  0  <  x  <  oo,  0  <  t 


u-v  =  g(t) , 

«  =  /"(*), 

v  =  d'(x) , 


i  =  0 


t  =  0. 


Corresponding  to  (2.6)  we  obtain  the  estimate 

ll*(*)HJ  +  \\m\\2  +  jf  (l“(°’r)|2  +  l»(0*T)|J)dr  <  C  (ll/l2  +  Ill'll2  +  jf  l5(r)|2dr)  . 

By  going  back  to  0  via  u,  v,  we  can  obtain  an  estimate  for  0,  but  now  in  terms  of  higher  order 
derivatives. 

As  we  demonstrated  above,  the  introduction  of  higher  derivatives  in  the  estimates  can  be 
avoided  by  integrating  the  boundary  condition  with  respect  to  time.  The  analogous  technique  can 
in  general  not  be  applied  for  difference  approximations.  Furthermore,  perturbations  in  derivative 
boundary  conditions  produce  a  possible  polynomial  growth  in  time,  which  propagates  into  the 
domain.  For  approximations  there  will  be  a  growth  of  the  solution  which  is  at  least  proportional 
to  n.  For  each  extra  factor  s  that  multiplies  the  function  B{u,s)  in  (2.2),  the  possible  polynomial 
growth  in  the  approximation  increases  one  order  in  n.  In  general,  one  should  therefore  avoid 
derivative  boundary  conditions,  for  first  order  systems. 

In  order  to  make  these  conjectures  more  precise,  we  consider  difference  approximations  to  the 
wave  equation.  Our  approach  is  to  use  a  technique  similar  to  the  one  for  the  continuous  case,  in 
order  to  obtain  estimates  of  the  solution.  We  treat  only  the  standard  5-point  formula,  and  consider 
first  the  Cauchy-problem 

■  -  2*?  +  <t>TX  =  A2(0?+ 1  -  20’  +  *7-1 ) .  A  =  k/h 

(2.11)  *7  =  fi 

where  k,  h  are  the  time-  and  space-steps  respectively.  (Other  initial  procedures  could  of  course  be 
used,  but  this  one  is  enough  for  our  purposes.)  After  Fourier-transformation  in  ar,  we  get 

*"+1  -  2  (1  -  2A2  sin2  ■£)  0"  +  0""1  =  0 , 

which  has  the  characteristic  roots  zj,  zj  on  the  unit  circle  if  A  <  1.  For  u;  =  0  they  coincide  at 
z  =  1;  accordingly  there  is  a  possible  growth  in  the  solution  of  order  n.  This  cannot  be  avoided 
with  any  consistent  method.  By  choosing  the  initial  data  such  that  u*  -  u®  =  0(k)  (corresponding 
to  a  bounded  d-function  in  (2.11)),  the  growth  is  at  most  of  order  tn  =  nk,  which  is  the  same  as 
for  the  differential  equation.  But  one  should  be  aware  of  the  fact  that  perturbations  in  the  initial 
data  are  magnified,  and  that  the  scalar  wave  equation  i  net  suitable  for  integration  over  very  long 
time  intervals. 

We  next  introduce  discrete  boundary  conditions.  We  shall  use  the  difference  operators 

A±x0"  =  =  ±(0"±1  -  0") 

A ±,*7  =  *z?±t0"  =  ±(0?±l  -  0") . 

The  inhomogeneous  version  of  the  first  condition  used  in  [1]  which  approximates  the  boundary 
condition  in  (2.5)  is 


(2.12) 


£+«(0O  +  0")  “  D+  r(0o  +  0Q  +  l)  =  2<7(tn+l/2)  • 


In  order  to  obtain  correct  centering,  the  grid  is  shifted  such  that  xo  =  -/i/2.  By  introducing  the 
new  variables 

<+i = > 

the  difference  scheme  can  be  written  in  the  form 


(2.13) 


The  boundary  condition  (2.12)  is 

(2.14)  ttf+1  +  t*5+1  -  (t£+1  +  tjJ )  =  2 g(tn+l/2) . 

From  now  on,  the  notation  ||  ■  ||  is  used  for  the  discrete  /2-norm  defined  by 


(a) 

u»+i  -  u*  =  A  , 

;  =  1,2,... 

(b) 

w»+i  -  =  A  A+,u"+1 , 

3  =  0, 1, . . . 

(c) 

ui  =  , 

II 

O 

(d) 

=  D+xfj , 

»— 4 

O 

II 

(2.15) 


n«i2 = i>?i  2fc- 


i= 0 


We  shall  prove 

Theorem  2.1.  The  boundary  condition  (2.14)  gives  a  stable  difference  scheme  (2.13)  for  A  <  1, 
and  for  0  <  t„  <  T  there  is  an  estimate 


(2.16)  ||«l2  +  ||t>"||2+£(Ki2 


Kl2)* 


£C( 


||d||24H|JWII2  +  I>(‘,-l/2)|2fc 


Proof.  The  normal  mode  analysis,  see  [6],  is  applied.  We  need  the  general  form  of  the  solution  to 
the  transformed  system 


(z-  1)6;  =  .7  =  1,2,... 

(z  -  1)6,  =  Az(ui+j  -  iij ) ,  ;'  =  0, 1 . 


When  requiring  bounded  solutions  for  |z|  >  1,  there  is  only  one  component  present,  and  it  has  the 
form 


(2.17) 


=  <7 


±K  y/z/*\ 


Here  the  constant  o  is  determined  by  the  boundary  condition,  and  k  satisfies  the  characteristic 
equation 


(2.18) 


(z-  l)2  =  A2  -(*-  l)2 
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W  z.j-.V.  c.Z.VlV.V.V.V.V.,-.  .  v 


•\  * 


i.e.,  we  have 


(2.19) 


-  1  =  ±A  y/i  («  - 


The  apper  and  lower  signs  in  (2.17)  and  (2.19)  correspond  to  each  other,  and  are  chosen  such  that 
|k|  <  1  for  |z|  >  1. 

The  transformed  boundary  condition  (2.14)  is 


*(«i  +  uo)  -  (z  +  1)  wo  =  2$ 


and  when  (2.17)  is  inserted,  we  get  for  g  =  0, 


(2.20) 


=  ±-J^(k 

k\z 


We  shall  prove  that  there  is  no  solution  z,  k  to  (2.19),  (2.20)  with  \z\  >  1.  By  squaring  both  sides 
of  (2.20)  we  get 

*2-(*+i)  *  +  i*o. 


Obviously  the  coefficient  of  *  in  this  equation  must  equal  the  coefficient  of  *  in  (2.18),  i.e., 


*+i  =  2(l-A2)  +  A2  (ie+ A)  , 


and  this  is  possible  for  A  <  1  only  if  k  ss  1.  This  corresponds  to  *  =  1  in  (2.19),  and  a  perturbation 
z  =  1  -M,  £  >  0  shows  that  the  corresponding  root  k  is  inside  the  unit  circle  if  the  lower  sign  is 
chosen.  But  in  that  case  (2.20)  does  not  hold,  showing  that  the  Kreiss  condition  is  fulfilled. 

The  norm  used  in  [6]  is  integrated  in  time,  furthermore  there  is  an  assumption  of  zero  initial 
data.  However,  if  the  approximation  fulfills  an  energy-estimate  for  periodic  solutions,  then  the 
/2-norm  of  the  solution  can  be  estimated  at  any  given  time  in  terms  of  the  initial  and  boundary 
data.  That  part  of  the  proof  is  omitted  here,  we  refer  to  [3]  for  the  proof  in  the  semi-discrete  case. 

The  energy  estimate  for  periodic  boundary  conditions  is  easily  obtained.  Define  the  scalar 
product  and  norm  for  scalar  periodic  grid-functions  by 


K  »)p  =  uivih '  iMIp  =  K  «)* 
1 


where  =  «j+/v.  By  taking  the  scalar  product  of  (2.13a)  with  iin+1,  and  of  (2.13b)  with  vn+1, 
wg  get 

ll«n+1|£  =  (un,ttn+1),  +  A(A_,rB,«"+1)p 

IK+1Hp  =  (vB,  vn+1)p  +  A(A+r«n+1,  un+l)p . 

From  (2.13a),  (2.13b)  we  also  get 


(un,un+1)p  =  ||u"||p  +  A(A_*t>n,  un),, 


(vn+\  Vn)p  =  ||„12  +  A(A+run+1,  wn)p  , 


and  by  using  the  identity 

(A_r»,u)p=  -(A+Xtt,  v)p 

we  obtain 

(2.21)  E(Un+1)  =  E(Un), 

where 

wn)  =  h%  +  \\v%  +  A(  A_,en,  u")p . 

Since 

|( A_ru,  u)p|  <  ||»||p  +  ||u||p , 

E(Un)  is  a  proper  energy-norm,  which  for  A  <  1  is  equivalent  to  ||un||2  +  ||vn||p.  This  proves  the 
theorem. 

By  using  the  definition  of  u  and  v  we  get  immediately  from  Theorem  2.1: 

Corollary  2.1.  The  solution  to  the  discrete  wave  equation  (2.11)  with  boundary  condition  (2.12) 
satisfies  the  estimate 

||*12  +  ||l?-t<n|2  +  \\D+x<t>n\\  +  l<Ao|2  +  I2  +  \D+Xc%\2)k 

K=1 

<  C  (W  +  \\D+Xf\\2  +  |/o|2  +  £  |s(f„-x/2)|2* 

V  11=1 

This  inequality  corresponds  exactly  to  the  corresponding  one  (2.8)  for  the  continuous  case. 

Next  we  shall  consider  the  higher  order  boundary  condition  (2.9)  for  the  one-dimensional 
scheme.  The  approximation  used  in  [1]  for  this  condition  is 

D+tD_t(<f>X  +*?)-*  D+x(<l>o+1  ~  *o-1)  =  M*n)  , 
which  stated  in  terms  of  the  variables  u,  v  takes  the  form 

(2.22)  t£+1  -  uS  +  u?+1  -  <  -  «+1  -  no""1)  =  2kg(tn) . 

The  stability  analysis  follows  the  same  lines  as  for  the  first  order  condition.  In  fact,  the  proof  of 
Theorem  2.1  is  identical  up  to  the  condition  (2.20),  which  for  (2.22)  becomes 

(2.23)  (z-l)[z+l*i«/f(«+l)J*0. 

k  y  is 

The  extra  factor  z  -  1  is  a  result  of  the  time-differentiation  applied  on  the  boundary  condition. 
Since  (2.23)  holds  for  z  =  1  regardless  of  the  sign  in  the  brackets,  there  is  necessarily  a  generalized 
eigenvalue  z  =  1  corresponding  to  it  =  1.  Therefore  the  scheme  is  unstable,  and  small  perturbations 
in  the  difference  scheme  or  the  boundary  conditions  can  be  expected  to  have  a  significant  effect  on 
the  solution. 


We  shall  further  investigate  the  approximation  in  order  to  see  if  it  is  possible  to  obtain  estimates 
under  more  restrictive  smoothness  assumptions.  This  will  explain  why  under  certain  circumstances 
higher  order  absorbing  boundary  conditions  can  be  used.  We  use  the  same  technique  as  indicated 
for  the  continuous  problem.  Corresponding  to  the  differentiated  variables  ut,  vt  we  define 

«”+1  =  D-tu]+1 ,  j  =  0,1,...;  n=  1,2,... 

S?+1  =  D-t  tf+1 ,  j  «  0, 1,  0, 1 . 

By  differencing  the  system  (2.13)  in  the  time-direction,  we  get 

|“i+1  =  AA-*5?,  j  =  1,2, ... ;  n  =  2, 3, . . . 

*■}  =  ■&(<% +<t>°j),  j  =  0,1,... 

=  &  (^j+1  -  *}  -  ^5+1  +  $)  >  j -0,1 . 

The  boundary  condition  (2.22)  takes  the  form 

uS+1  +  i?+1  ~  (®?+1  +  ®o )  =  2 g(tn)  ■ 

Apparently  un,  vn  satisfy  the  same  difference  equation  and  the  same  type  of  boundary  conditions 
as  the  solution  (un,un)  to  (2.13),  (2.14)  does.  Therefore,  an  estimate  of  the  type  (  2.16)  holds. 
However,  in  this  case  the  initial  data  it2,  t;1  as  given  in  (2.24)  are  different,  and  this  is  an  important 
distinction.  Even  if  the  grid-function  <f>  is  bounded,  it  cannot  be  expected  to  be  smooth  unless  we 
are  very  careful  when  initializing  the  problem.  In  particular,  it  may  be  difficult  to  match  the  initial 
data  to  the  boundary  data.  If  this  is  not  done  properly,  accuracy  will  be  lost.  (In  [1]  the  need  for 
compatible  boundary  and  initial  data  is  pointed  out,  but  the  connection  to  stability  and  accuracy 
is  not  discussed.) 

As  an  example,  consider  a  case  where  <£(x,  t)  =  0  for  t  <  0,  and  where  a  boundary  function  g(t) 
begins  driving  the  solution  at  t  =  0,  i.e.,  g(t)  =  0  for  t  <  0.  We  compare  the  boundary  condition 

(2.25)  <f>t-<t>,=9(t),  z  =  0 

with  its  differentiated  version 

4>tt  -  <t>xt  =  $'(<)  .  2  =  0, 

corresponding  to  the  first  and  second  order  conditions  (2.4)  and  (2.9)  respectively.  For  convenience 
we  extend  the  grid  one  step  backwards  in  time,  such  that  the  initial  conditions  are 

*-‘=*9=0,  j  =  0, 1 . 

The  discrete  first  order  boundary  condition  is  (2.12).  We  use  the  first  order  system  (2.13)  to 
compute  the  first  points.  It  turns  out  that  A  =  1  gives  the  best  local  accuracy  near  x  =  0,  t  =  0, 
and  a  direct  calculation  shows  that  with  (2.13),  (2.14),  we  get 

u\=g(k/2)  =  g(0)+$g'(0)  +  O(k1). 


Since  the  left-going  wave  u  +  v  is  zero,  the  true  solution  is  u(x,  t)  =  Q.5g(  t  -  x)  for  x  <  t.  Recalling 
that  ti"  is  centered  at  (xj,  fn— 1/2 )» the  true  solution  corresponding  to  «i  is 

u(k/2,3k/2)  =  0.5p(Jfc)  =  0.5[5(0)  +  kg'(O)  +  0(k2)} . 

Apparently  a  second  order  error  is  obtained  only  if 

9(0)  =  0. 


The  conclusion  is  that  if  g'( 0)  ^  0,  the  approximation  becomes  worse  with  the  higher  order  bound¬ 
ary  condition.  Both  boundary  conditions  are  centered  properly,  so  this  is  an  illustration  of  the  fact 
that  higher  order  derivatives  in  the  boundary  conditions  require  smoother  solutions.  Convergence 
may  still  occur  as  h  -*  0,  which  was  demonstrated  for  a  scalar  equation  in  [5],  but  in  general 
accuracy  is  lost. 

A  simple  numerical  experiment  was  performed  in  order  to  illustrate  the  results  obtained  above. 
We  used  the  interval  0  <  x  <  1,  and  the  right  boundary  was  treated  similarly  to  the  left  one,  but 
with  zero  data.  With  zero  initial  data  /(x),  d(x),  the  solution  is  generated  by  the  boundary  data 
g(t)  in  (2.5),  and  the  true  solution  d>(x,  t)  is  easily  derived. 

The  three  alternative  boundary  conditions  are 


BCO 

BCl 

BC2 


(W  =  4>(-h/2,tn) 

U&+,  =  0<l  +  /»/2,tn) 

I  A +t(<fiS  +  4>?)  -  A  A +x(<fiS  +  <t>S+l)  =  2 kg(tn+t/2) 

\  A+t(<fiy  +  <f>nN+l)  +  A  A +x(<t>nN  +  4>Vl )  =  0 

f  A+tA-t(4>Q  +<£")-  A  A+X(0o+1  -  4>o-1 )  =  2kg'{tn) 
\  A+tA )  +  A  A+x(^+1  -  <t>nN~l )  =  0  . 


For  each  boundary  condition  three  cases  were  run: 

NOPER:  No  perturbation 

INPER:  Perturbation  0.001  in  ,  j  =  0, 1, . . . ,  N  +  1 

BOPER:  Perturbation  0.001  in  the  right  boundary  condition. 

In  all  cases  the  boundary  function  g(t)  =  sin  t  was  used,  which  gives  a  discontinuity  in  the 
second  derivatives  of  4>.  Table  2.1  shows  the  mean-square  error  in  the  solution  at  t  =  4. 


NOPER 

INPER 

BOPER 


B< 

h  =  1/20 

:o 

h  -  1/40 

B< 

h  =  1/20 

:i 

h  =  1/40 

B< 

h  =  1/20 

22 

h  =  1/40 

7.0e  -  5 

3.8e  -  3 

6.0e  -  4 

2.1e  —  5 

3.9e  -  3 

4.2e  -  4 

2.0e  -  4 

1.4e  -  2 

4.0e  -  2 

5.1e  -  5 

2.6e  -  2 

7.9e  -  2 

3.5e  -  2 

1.4e  -  1 

3.5e  +  0 

1.8e  -  2 

2.2e  -  1 

1.4e  +  1 

Table  2.1.  Mean-square  error  in  <f>  at  t  =  4. 


The  perturbations  have  a  disastrous  effect  for  BC2,  in  particular  when  the  perturbation  occurs  in 
the  boundary  condition.  BC1  is  also  more  sensitive  to  perturbations  than  BCO,  as  expected.  When 
no  perturbation  is  introduced,  BC2  still  gives  a  poor  accuracy,  and  the  convergence  rate  is  clearly 
only  of  first  order. 

We  have  not  included  any  forcing  function  in  the  discussion  above.  When  deriving  error- 
estimates  for  the  approximation  0",  such  a  forcing  function  F"  represents  the  truncation  error,  and 
is  necessarily  present  in  the  difference  scheme  for  the  error  <p"  -  <£(ij,  tn).  The  grid  function  F" 
depends  on  </>(x,  t),  and  is  small  only  if  <$>  is  sufficiently  smooth.  The  results  above  show  that  if  this 
is  not  the  case,  the  boundary  condition  BCl  is  more  forgiving  than  BC2. 

For  multi-dimensional  problems  it  is  more  difficult  to  keep  the  solution  smooth  enough,  in  par¬ 
ticular  when  the  numerical  method  requires  extra  numerical  boundary  conditions.  The  conclusion 
from  this  section  is  therefore  that  derivative  bundary  conditions  should  be  avoided  for  hyperbolic 
problems.  In  the  next  section  we  shall  show  how  accurate  and  well-posed  conditions  still  can  be 
derived. 

3.  First  order  systems 

In  this  section  we  consider  general  first  order  hyperbolic  systems.  We  shall  discuss  the  im¬ 
plementation  of  the  boundary  conditions  and  how  to  provide  accurate  data  at  the  boundary.  If 
the  system  has  variable  coefficients  or  is  non-linear,  then  the  coefficients  are  frozen  outside  the 
computational  domain;  usually  the  values  at  infinity  are  used.  It  is  important  to  distinguish  be¬ 
tween  this  procedure,  and  the  method  of  using  values  from  the  state  at  infinity  for  the  variables 
at  the  computational  boundary.  In  the  first  case,  the  solution  is  still  allowed  to  vary  outside  the 
computational  domain,  but  it  satisfies  a  simpler  differential  equation.  In  this  way  the  form  of  the 
solution  can  be  explicitly  derived,  such  that  boundary  conditions  can  be  obtained.  As  in  Section  2. 
we  shall  consider  the  case  where  the  computational  domain  is  infinite  in  the  y-direction  but  has  a 
boundary  at  x  =  0.  The  right  half-plane  is  the  domain  for  computation,  the  solution  in  the  left 
half-plane  is  supposed  to  have  a  simpler  structure,  such  that  this  part  can  be  eliminated  from  the 
computation.  (In  practice,  there  is  of  course  another  computational  boundary  for  some  positive  x, 
as  well  as  in  both  y-directions.)  In  order  to  derive  conditions  at  x  =  0,  we  therefore  consider  the 


problem 

((a) 

Ut  +  AUr  +  BUy  =  0 , 

-00  < 

x  <  0,  -oo  <  y  <  oo,  0  <  t 

(3.1) 

{  (b) 

suPr,y  W(x,y,t)\  <  00  f 

0  <  t 

1(c) 

U{x,y,  0)  =  f(x,y) , 

-00  < 

x  <  0,  -oo  <  y  <  oo  . 
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In  [4]  we  have  treated  this  problem,  and  we  briefly  summarize  the  results  here.  It  is  a  generalization 
of  the  absorbing  boundary  conditions  derived  in  [1]. 

Since  the  system  is  hyperbolic,  the  matrix  A  can  be  diagonalized,  and  it  is  assumed  that  this 
has  already  been  done  here,  such  that 

A  =  diag(ai,a2,. .  .,a(+m) 

where 

>  0 ,  j  —  1, 2,  • . . ,  / 
aj+i  ^  Oi  j  —  1,2, 

The  system  (3.1a)  is  Laplace-transformed  in  time  and  Fourier-transformed  in  the  y-direction. 
s  and  w  being  the  dual  variables.  For  the  transformed  vector  U  we  get  the  equation 

=  -sQ{£)U(x,u,s)  + A~l  f{x,uj). 


where 

Q(0  =  A-l(/  +  {fl),  4  =  tw/e, 

and  where  /(x,w)  is  the  Fourier-transformed  initial  function.  For  Re  s  >  0,  the  matrix  0(4)  can 
be  transformed  to  block-diagonal  form 


T-l(tX)(  or(0* 


Q*(t) 

o 


o 

Q"(0  ’ 


where,  corresponding  to  the  eigenvalues  of  A,  07(4),  QlI(0  are  (/  x  /)-  and  (m  x  m)-matrices 
respectively.  The  condition  (3.1b)  implies  the  condition 


(3.2) 


(T-^O^O.^s)]'  =  -  jf°  f(a,u)]1  da 


where  [V]*  indicates  the  upper  /  elements  of  the  vector  V.  This  is  the  exact  boundary  condition 
for  the  problem  defined  in  the  domain  i  >  0.  For  /  s  0,  it  is  the  completely  absorbing  boundary 
condition  by  Engquist  and  Majda.  In  order  to  get  local  boundary  conditions  in  physical  space, 
some  approximations  must  be  made. 

Let  Kj,  j  =  1,2, . . .,/  be  the  eigenvalues  of  -s07(4)i  and  assume  that  Ql(0  can  be  transformed 
to  diagonal  form.  Assuming  that  this  transformation  is  included  in  the  matrix  T(4)  above,  the 
equation  (3.2)  can  be  written  in  the  form 

(3.3)  (T-1(O]>&r(0,u;,s)=-(T-1(OA-1]i  f°°  e~H^  f (a, u)d<T ,  j  =  1,2 . / 

Jo 

where  [A]j  denotes  the  j-th  row  of  the  matrix  A.  If  |4|  is  small,  we  use  the  expansions 


(3.4) 


(a)  r-*(4)  =  /+*4  +  54*  +  0(|4l3) 

(b)  a,  =  -s(a-‘  +  0#  +  0(|4|J)) ,  J  =  1.2- . / 
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in  (3.3).  If  £- terms  of  order  one  and  higher  are  disregarded,  we  get  after  transforming  back  to 
physical  space 

(3.5)  U^(0,y,t)  =  ^(~ajt,y),  j  =  l,2,...,/. 

This  condition  is  obtained  directly  from  (3.1a)  by  letting  B  =  0,  and  then  solving  the  one¬ 
dimensional  system  exactly.  If  the  £-terms  of  first  order  are  kept,  we  get  after  a  straightforward 
calculation 

(3.6)  lf(i)( 0,  y,  t)  +  /  RM 0,  y,  t)  dr  =  g(y,  t) , 

Jo 

where 

ff(y>0  =  +  /  fy(—ajT,y  ~  0jajT)  dr ,  j  =  1,2,...,/. 

Jo 

These  are  the  boundary  conditions  for  the  problem 

{Ut  +  AUX  -I-  BUy  =  0  ,  0  <  x  <  oo,  -oo  <  y  <  oo,  0  <  t 

saPr,y  <  °°,  0<* 

U ( x ,  y,  0)  =  /(z,  y),  0  <  x  <  oo,  -oo  <  y  <  co. 

There  is  no  proof  of  well-posedness  for  general  systems,  but  we  assume  that  (3.7),  (3.6)  is  well-posed 
such  that  there  is  an  estimate 

(3-8)  ||tf (•, ., t)||2  +  j f*  \\U(0,  -,r)||2  dr  <c  (||/(-,  -)l|2  +  f  ||*(-,  0||2  dr)  . 

Here  ||f/(*, -,/)||  is  the  Z2-norm  over  the  right  half-plane,  and  ||f/(x,  -,t)||  is  the  J^-norm  over 
{y/  -  oo  <  y  <  oo}. 

At  this  point  we  make  a  remark  concerning  the  theory  for  well-posedness.  The  Kreiss  theory 
was  developed  for  first  order  systems  with  boundary  conditions  of  the  type 

SoU{0,  y,t)  =  y(y,t) 

where  So  is  a  constant  rectangular  matrix.  In  transformed  space  the  boundary  condition  can  be 
written  as 

H  V,*)V'(0,u;,3)  =  B"{u,  s)Vn(0,u,  s)  +  y(w, s) ,  V  =  T~l  U , 

and  the  condition  for  well-posedness  is 

(3.9)  Det(jB/(w,s))  ^  0 ,  Imu/  =  0,  Res>0. 

For  s  ^  0,  it  is  clear  from  the  construction  above,  that  B1  is  a  function  of  the  single  variable 
£  =  iu/s.  For  uj  —  0,  B1  is  a  constant  matrix,  consisting  of  the  first  /  columns  of  So,  accordingly, 
the  condition  for  well-posedness  is  independent  of  s.  It  follows  from  these  arguments  that  it  is 
sufficient  to  consider  the  domain  |w|  +  |s|  =  1  when  verifying  the  Kreiss  condition. 

If,  on  the  other  hand,  derivatives  occur  in  the  boundary  conditions,  then  B1  may  depend  on  s 
even  for  u>  =  0,  and  the  point  u  =  s  =  0  must  be  included. 

If  the  conditions  (3.6)  are  differentiated  with  respect  to  t,  we  get 

(3-10)  V}j)(0,y,t)+  RjUv(0,y,t)  =  gt(y,t),  j  =  1,2 . /. 

These  conditions  are  less  suitable  for  computations.  We  have 
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Lemma  3.1.  The  problem  (3.7),  (3.10)  is  at  least  weakly  ill-posed. 

Proof.  For  u>  =  0,  the  condition  (3.9)  is 

Det(s/)  ^  0  ,  Re  s  >  0 
which  is  obviously  violated  at  s  =  0. 

As  for  the  wave  equation,  we  conclude  that  the  condition  (3.6)  is  a  better  basis  for  an  ap¬ 
proximation  than  the  condition  (3.10).  In  [4]  an  experiment  was  carried  out,  clearly  demonstrating 
this.  Note  that  the  condition  (3.6)  is  no  more  difficult  to  implement  than  the  condition  (3.10).  The 
integral  is  substituted  by  a  sum,  and  for  each  y- value  one  new  term  is  added  for  each  times-step. 

In  many  applications  there  is  a  source  outside  the  computational  domain,  which  is  driving  the 
solution.  Sometimes  this  situation  can  be  modeled  as  an  initial  value  problem  where  the  initial 
function  /(x,  y)  is  known  everywhere.  However,  in  order  to  avoid  the  approximations  that  are 
necessary  in  the  general  case,  a  better  method  may  be  to  prescribe  the  correct  data  directly  at  the 
boundary. 

Assume  that  the  solution  consists  of  two  parts  Uq,  U\ ,  where  Uq  is  known.  A  typical  case  is 
where  Uq  is  a  plane  wave  produced  by  an  outer  source.  In  the  computational  domain  there  is  some 
interaction  which  causes  the  generation  of  U\,  and  it  is  assumed  that  this  part  of  the  solution  can 
be  considered  as  a  wave  that  passes  out  through  the  boundary.  In  such  a  case,  the  initial  data 
occurring  in  (3.5)  and  (3.6)  are  substituted  by  the  known  forced  solution  Uq.  The  conditions  (3.5) 
and  (3.6)  for  U  -  Uq  +  U\  takes  the  form 

(3.11)  U<»(0,  y,  t)  =  U^(0,  y,t),  >  =  1,2 . / 


(3.12)  £f(j)(0,y,t)+  f  RjUv(0,y,T)dT=U(oi)(0,y,t)+  fjl;[(tro),](0,y,r)dr,  j  =  1,2,...,/. 
Jo  Jo 

Consider  next  the  wave  equation  (2.1).  From  the  discussion  in  Section  2  it  is  dean*  that  there  are 
certain  advantages  with  using  the  corresponding  first  order  system  as  a  basis  for  the  approximation. 
With 
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(3.13) 

vt  + 

1 

0 

0 

V , 

:  + 

0  0 
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Vv  =  0. 


The  coeffident  matrix  multiplying  Vx  is  singular,  hence  the  technique  used  above  does  not  apply 
directly.  We  modify  the  system  by  subtracting  a  term  eKr,  (  >  0,  and  after  deriving  the  boundary 
conditions  we  let  c  — ►  0.  With  the  new  variables 

u  +  v 


u  = 


V2  ’ 


v  — 


U  -  V 

7T’ 


w  =  w , 
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we  get  the  system 

(3.14) 

where 


Ut  +  AUX  +  BUy  =  0 


A  =  0  — 1  — «  0  ,  B  =  —=  0  0  1- 

V2 

L  o  o  -d  Li  i  oJ 


0  0  n 


The  matrix  R  and  the  eigenvalues  Kj  in  (3.4)  are 


R=7i  0  0  -  1 

Li  - 1  i  +  f  o . 

=  ~3  (rb  ~  V  +  °(l^l4)) 

«2  =  -*  (“T+7  +  V  +  ^(l^l4)) 

«3  =  j/<  . 

There  is  only  one  eigenvalue  kj  with  Re  k  <  0  for  Re  a  >  0,  hence  the  boundary  condition  is 
determined  by  the  first  row  of  R.  The  condition  (3.12)  becomes 

(3.15)  r(0,y,t)  +  wv(0,y,T)dr  =  «(-(l  -  c)t,y,0)  -  J  tnv(-(l  -  f)r,y,0)dr. 

Since  the  matrix  row  in  (3.6)  is  well  defined  even  for  e  =  0,  the  limit  of  (3.15)  as  e  —  0  is 

also  well  defined.  In  the  original  variables  the  condition  becomes 


(3.16) 


(**  +  ®1(0,  y,  t)  =  [u  +  »](  — t,  y, 0)  —  /  wv(-r,y,0)dr. 

Jo 


If  the  right  hand  side  is  considered  as  given  data,  the  condition  (3.16)  is  trivially  well-posed.  We 
note  that  for  zero  initial  data,  there  is  no  difference  between  the  first  and  second  order  conditions. 
A  third  order  condition  would  be  obtained  by  including  the  third  term  in  the  expansion  (3.4a), 


+  5a- 


l-(2 


’1  1  O' 

110, 

.0  0  0. 


Considering  zero  initial  data,  we  get  from  the  transformed  differential  equation  (3.14) 


£(«  +  v)  =  \/2  wx  -  tfc)  , 


and  by  inserting  this  into  the  second  row  of  5,  we  get  for  «  =  0, 


(3.17) 


u  +  -i*w  =  0. 

2v/2 


In  physical  space  this  condition  is  for  the  original  variables 

1  f* 

(3.18)  [«  +  »](0,y,t)+ -  J  «j,(0 ,y,r)dr  =  g(y,t), 

where  we  have  included  the  possibility  of  giving  non-zero  boundary  data.  By  differentiating  with 
respect  to  t  and  going  back  to  the  wave-function  <t>  in  (2.1),  we  obtain 

(3.19)  —fat  +  <t>xt  +  ^ <t>n  =  0t(lM)  *  *  =  0, 

which  we  recognize  as  the  condition  (2.9)  with  g  =  —  gt. 

The  well-posedness  of  the  first  and  second  order  conditions  is  trivial.  For  the  third  order  one 
we  have 

Theorem  3.1.  The  wave  equation  formulated  as  a  first  order  system  (3.13)  is  well-posed  with  the 
boundary  condition  (3.18). 

Proof.  The  first  column  of  the  matrix  T  which  transforms  Q{£)  to  diagonal  form  is 
[T(0](l)  =  [l-*i/a,  1  +  Kt/a,  — >/2(l  -  k?/'*2)/£]T , 

Referring  to  Section  4  of  [4],  we  obtain  the  condition  for  well-posedness  from  (3.17)  as 


*t  =  ±e  yjl  -  &  . 


The  sign  of  *i  is  taken  such  that  Re  <  0  for  Re  s  >  0.  (Note  that  it  is  not  sufficient  to  consider 
small  values  of  |£|  here.)  The  critical  points  are  therefore  given  by 

(3.20)  1  T  a/1  =  0 

and 

(3.21)  1-^  (l±  Vl-f2)  =0. 

The  condition  (3.20)  implies  (  =  0,  which  is  the  one- dimensional  case  known  to  be  well- posed.  The 
condition  (3.21)  implies 

±vi-e  =  i, 

which  again  leads  to  the  one-dimensional  case  £  =  0,  and  the  theorem  is  proven. 

If  one  prefers  to  use  the  scalar  wave  equation,  then  the  condition  (3.18)  takes  the  form 

-^«(0,  y, t)  +  ^*(0, ¥,<)+$■  /  fa(0,y,r)dr  =  g(y,t) 

Jo 


which  is  the  integrated  form  of  (3.19).  In  one  space-dimension  it  is  equivalent  to  the  first  order 
condition  (2.4)  denoted  by  BCl  in  Section  2;  the  condition  (3.19)  corresponds  to  BC2.  Table  2.1 


shows  that  the  integrated  version  is  preferable.  The  numerical  experiments  prescribed  in  [4]  for 
two-dimensional  problems  confirms  this  conclusion. 
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